On coarse- graining by the Inverse Monte Carlo method: Dissipative Particle Dynamics simulations 

made to a precise tool in soft matter modeling 
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We present a promising coarse- graining strategy for linking micro- and mesoscales of soft matter systems. 
The approach is based on effective pairwise interaction potentials obtained from detailed atomistic molecular 
dynamics (MD) simulations, which are then used in coarse-grained dissipative particle dynamics (DPD) simu- 
lations. Here, the effective potentials were obtained by applying the Inverse Monte Carlo method [Lyubartsev 
& Laaksonen, Phys. Rev. E. 52, 3730 (1995)] on a chosen subset of degrees of freedom described in terms of 
radial distribution functions. In our first application of the method, the effective potentials were used in DPD 
simulations of aqueous NaCl solutions. With the same computational effort we were able to simulate systems 
of one order of magnitude larger as compared to the MD simulations. The results from the MD and DPD 
simulations are found to be in excellent agreement. 
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I. INTRODUCTION 

Classical molecular dynamics computer simulations with 
site-site pair potentials can currently be run for molecular sys- 
tems consisting of the order of 10 4 atoms, corresponding typi- 
cally to a system size of 50-60 A. Corresponding time scales 
using the same approach are typically of the order of 10 ns. 
This may be enough for simulations of isotropic liquids of 
simple molecules, while for many complex systems the re- 
quirements for the length and time scales are much more ex- 
tensive. 

The restrictions of the detailed all-atom description are es- 
pecially severe in soft matter systems such as biomembranes 
and polyelectrolytes, although these are just two examples of 
the diversity of soft materials. This is mainly due to a wide 
range of time and length scales associated with these systems, 
and it is obvious that these problems become particularly pro- 
nounced in studies of dynamic processes which usually take 
place under hydrodynamic conditions. The primary aim in 
coarse-graining is to build a bridge between different time 
and length scales. Although there is no unique way to do 
this, the coarse-graining procedure always leads to a reduction 
of the number of degrees of freedom and thus alleviates the 
computational burden. The challenge is to establish system- 
atic coarse-graining schemes that allow one to develop sim- 
plified model systems in terms of the information extracted 
from the underlying microscopic systems. Thus far, various 
approaches have been used to coarse-grain atomic and molec- 
ular systems. Below, we will briefly discuss some of them. 

As a relatively generic approach, the dissipative particle dy- 
namics (DPD) | Jyl HI has been used to study coarse-grained 
particle and polymeric systems, see e.g. Ref. |4] for a review. 
In our approach we use DPD as a thermostating method which 
is discussed in more detail in the following sections. To men- 
tion other relatively generic approaches, Flekk0y et al. 



have recently introduced a hybrid method to combine particle 
and continuum level descriptions. This framework allows a 
derivation of the DPD model and in practice uses a Voronoi 
tessellation based technique for simulations of coarse-grained 
particles at different levels. Another very promising approach 
was presented by Malevanets and Kapral QSl wno proposed 
a framework to couple a molecular level description with a 
mesoscale treatment of solvent that conserves hydrodynam- 
ics. This approach seems particularly appealing for studying 
dilute systems with hydrodynamics. 

A lot of impetus for multiscale modeling has come from 
polymer research due to both fundamental and practical rea- 
sons. Akkermans and Briels |2I applied a projector operator 
formalism based approach to coarse-grain polymer chains us- 
ing a microscopic model as a starting point. Their approach 
is appealing and deserves further studies since it accounts for 
the fluctuating forces at the coarse-grained level. Bolhuis et 
al. fiot Hill used the Ornstein-Zernike equation fl2ll with a 
hypernetted-chain closure to invert the pair correlation func- 
tion in order to derive effective pair potentials for neutral poly- 
mers in a solution. The coarse-grained polymers were repre- 
sented as soft colloidal particles. Their approach is based on 
a theorem proved by Henderson in 1974 II 311 that stipulates 
a unique, point-to-point correspondence between pair poten- 
tials and radial distribution functions g(r). The importance of 
this theorem becomes clear when one notices that for any g(r) 
under given conditions there is a unique pair potential that in- 
cludes corrections from the many-body interactions, and that 
by definition the pair potential then yields the original g(r). 
The problem is thus reduced to inverting the g(r) obtained 
from a microscopic model - how to obtain a coarse-grained 
model from that information lies at the heart of the problem 
of coarse-graining. At the representational level the soft el- 
lipsoid model for polymer melts by Murat and Kremer itbUl 
is very close to the one of Bolhuis et al. fiol Hill , but in- 
stead of using the pair correlation function Murat and Kremer 
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base their approach on separating the free energy into intra- 
and interchain parts. In an attempt to model more realistic 
systems, Reith et al. |15] introduced a coarse-graining pro- 
cedure that starts from a detailed chemical description of a 
polymer. In order to obtain a coarse-grained approach, a sim- 
plex algorithm based optimization is performed to obtain a 
coarse-grained force field. The structural features produced 
by this approach are in good agreement with experimentally 
observed ones. For a more complete overview of the current 
coarse-graining approaches, see e.g. Ref. 11611 and references 
therein. 

Our approach is also based on inverting the radial distri- 
bution function but instead of the integral equation based ap- 
proach of Bolhuis et al. we use the inverse Monte 
Carlo procedure of Lyubartsev and Laaksonen 1 17] which will 
be discussed in the following section. This systematic coarse- 
graining method is then combined with a momentum conserv- 
ing dissipative particle dynamics (DPD) thermostat. Our ap- 
proach consists of three conceptually simple steps. First, we 
use a detailed atomistic description and perform MD simu- 
lations. From these simulations we compute radial distribu- 
tion functions g(r) between different atoms, molecules, or 
molecular groups. Second, we apply the Inverse Monte Carlo 
ifTvtl procedure to invert the radial distribution functions. This 
yields effective interaction potentials V (r) between the se- 
lected interaction sites. Finally, we employ these effective, 
coarse-grained potentials in the DPD algorithm to study the 
large-scale properties of systems in question. Note that since 
DPD satisfies momentum conservation, the present approach 
allows studies of a given system with full hydrodynamics. 

Preliminary results of this work were recently reported in 
Ref. d, 

in which we presented only the most essential ideas 
of this approach. Our aim in this article is to provide one 
with a comprehensive discussion of issues related to coarse- 
graining and mesoscopic simulations in terms of the DPD 
method. In particular, we discuss in detail how atomistic 
MD simulations can be coupled to DPD by the Inverse Monte 
Carlo method. To this end, we use NaCl as our model system. 
Even in this simple system, and with a very modest degree of 
coarse-graining, we obtain a computational speed-up of one 
order of magnitude. Our procedure is well-defined and truly 
allows easy and controllable tuning of the level of coarse- 
graining. Finally, we show through coarse-grained DPD sim- 
ulations that this approach is physically sound and provides 
results in excellent agreement with the MD simulations and 
experiments. 

This paper is organized as follows. First, in Sect. |n] we 
describe the methods: The Inverse Monte Carlo procedure and 
the DPD algorithm. In Sect. |in| we present the results and 
finally, in Sect. II VI we discuss the results and some general 
aspects related to the general applicability of the method. 



II. METHODS AND MODELS 

As discussed in the Introduction, detailed atomistic simu- 
lation techniques are suitable for studies of microscopic fea- 
tures of soft matter systems, e.g. the excited states of a chro- 



mophore where the time scales of interest are within a few 
hundreds of femtoseconds. For studies of systems that are 
characterized by large time and length scales, such as a pro- 
tein in water solution over tens or hundreds of nanoseconds, 
however, atomistic simulation techniques are not very feasi- 
ble. Our aim in this section is to discuss ways how this prob- 
lem can be resolved by coarse-graining microscopic systems. 
In particular, we concentrate on an approach based on solv- 
ing the "inverse problem" which yields effective interaction 
potentials which in turn can be coupled to mesoscopic sim- 
ulation techniques such as the dissipative particle dynamics 
method discussed below. Full details of the methods can be 
found in Refs. 1 17] and 1 19, 20], respectively. 



A. The Inverse Problem in statistical mechanics 

In the standard statistical mechanical description of soft 
matter systems one begins from a formulation of a model, 
which is usually given in terms of a Hamiltonian or interac- 
tion potentials between atoms and molecules. When a Hamil- 
tonian is specified, one can use different numerical methods 
to compute canonical averages such as energies, distribution 
functions, and response functions which can be compared 
with experiments. As one has an explicit expression for the 
canonical averages, their calculation is in principle a fairly 
straightforward task, although in many cases it may be com- 
putationally a very demanding one. 

Occasionally one is interested in the solution of the inverse 
problem, i.e., how to deduce information about molecular in- 
teractions from the canonical averages which may be known 
from experiments. More specifically, the question is how to 
reconstruct the Hamiltonian based on results for the canon- 
ical averages. Many experimental properties can be chosen 
as a starting point to this end. Among the most important 
ones for the liquid state are the radial distribution functions 
g(r), which may be known from the structure factor measure- 
ments using X-ray or neutron scattering |T2j|. In 1974, Hen- 
derson proved a theorem [13] stipulating a unique, point-to- 
point correspondence between pair potentials and radial dis- 
tribution functions. In essence it states that for a given system 
under given conditions of temperature and density, two pair 
potentials which give rise to the same radial distribution func- 
tion cannot differ more than by an additive constant. As this 
constant can be defined from the condition that the interaction 
potential tends to zero at an infinite distance, the potential is 
uniquely defined. In practice the solution of an inverse prob- 
lem is not a straightforward calculation, however, since g{r) 
does not provide one with a formal expression for the poten- 
tial. Therefore, some special techniques are required for the 
solution of the inverse problem. 

In 1988, McGreevy and Pusztai suggested a Reverse Monte 
Carlo (RMC) method in which the starting point for a simula- 
tion was the radial distribution function |21]. In this method, 
Monte Carlo simulations are carried out without any prior 
knowledge of the interaction potential to fit g(r) that serves 
as an input. The set of configurations obtained in RMC may 
be used for further structural analysis, for example for calcula- 
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tion of three-dimensional spatial or orientational correlations. 
However, as the interaction potential is not recovered, the in- 
verse problem is not solved completely, and it is not possible 
to calculate thermodynamical or dynamical properties using 
this approach. For the same reason this approach in not suit- 
able for the development of coarse-grained models. 

In another approach, Reatto et al. 1 22] used the Hypernet- 
ted Chain (HNC) approximation to solve the inverse problem. 
This and similar approaches, based on some closure of the 
Ornstein-Zernike equation, have been used in a number of 
works during the last decade to compute interaction potentials 
from radial distribution functions 1 231 12411 . It should be noted, 
however, that computations within the HNC theory are fea- 
sible only for relatively simple models. Moreover, although 
yielding sometimes very accurate results 1251 . the HNC the- 
ory is not an exact mathematical solution of the statistical- 
mechanical problem, and its accuracy should be investigated 
carefully in every specific case. Other works devoted to the 
inverse problem are listed in Refs. I! 2al27l 12811 . 

A practical way to solve the inverse problem has been sug- 
gested by Lyubartsev and Laaksonen flm . This is the method 
we will concentrate on in the following, and we will use it to 
compute effective potentials from radial distribution functions 
obtained from detailed MD simulations. The effective poten- 
tials then allow us to to study the large-scale properties of a 
given model system through coarse-grained DPD simulations. 



B. The Inverse Monte Carlo Method - A Simple Case 

To illustrate the Inverse Monte Carlo method, we consider 
the case of a single-component system consisting of identi- 
cal particles with pairwise interactions. The general case of 
a multicomponent system is a straightforward extension of it. 
The Hamiltonian for this system is given as 



(i) 



where V(rij) is the pair potential and r,j is the distance be- 
tween the interaction sites i and j. Let us assume that we 
know the radial distribution function g(rij). Our aim is now 
to find the corresponding interaction potential Vfaj). 

Let us introduce the following grid approximation to the 
Hamiltonian, 



V(r) = V(r a ) = V a 



(2) 



for 



(3) 

where a = 1, . . . , M, and M is the number of grid points 
within the interval [0,r cut ], and r cut is a chosen cut-off dis- 
tance. Then, the Hamiltonian in Eq. Q can be rewritten as 



H 



(4) 



where S a is the number of pairs with interparticle distances 
inside the a-slice. Evidently, S a is an estimator of the radial 
distribution function: (S a ) = Anr 2 : g(r a )N 2 / \2V). The av- 
erage values of S a are some functions of the potential V a and 
can be written as an expansion 



(5) 



The derivatives d(S a ) / dV 1 can be expressed using exact sta- 
tistical mechanics relationships 1 17] 



d(S a ) d /dq5 a ( g )exp(-/3EA^A^A(g)) 

dV 1 8V 1 Jdqexp(-/3J2 x K x S x (q)) 



(6) 



Equations (|5j and allow us to find the interaction poten- 
tial V a iteratively from the radial distribution functions (S a ). 
Let Va ^ be a trial potential for which a natural choice is the 
potential of mean force (PMF) 



V<® = -k B Tlng(r a ), 



(7) 



By carrying out standard Monte Carlo simulations, one can 
evaluate the averages (S a ) and their deviations from the ref- 
erence values S 1 * defined from the given radial distribution 
function as A{S a )^ = (S a )M - S*. By solving the sys- 
tem of linear equations [Eq. Q] with coefficients defined by 
Eq. (|6ji, and omitting terms 0(AV 2 ), we obtain corrections 

to the potential A vi°^ . The procedure is then repeated with 

the new potential Va = vi°^ + AV'i ' 1 until convergence is 
reached. The whole procedure resembles a solution of a mul- 
tidimensional non-linear equation using the Newton-Raphson 
method. 

If the initial approximation of the potential is poor, some 
regularization of the iteration procedure is needed. In this 
case we multiply the required change of the radial distribu- 
tion function by a small factor that is typically between and 
1. By doing so, the term 0(AV 2 ) in Eq. Q can be made 
small enough to guarantee convergence of the whole proce- 
dure, although the number of iterations will increase. 

The above algorithm has also the advantage that it provides 
us with a method to evaluate the uncertainty of the inverse 
procedure. The radial distribution function has normally some 
uncertainty. An analysis of the eigenvalues and eigenvectors 
of the matrix in Eq. allows one to make conclusions of 
which changes in g(r) correspond to which changes in the 
potential. For example, eigenvectors with eigenvalues close 
to zero correspond to changes in the potential which have al- 
most negligible effect on the radial distribution function. The 
presence of these small eigenvalues makes the inverse prob- 
lem not well-defined, however. In some cases, such as for 
liquid water, that may pose serious problems in the inversion 
procedure 12911 . 

It is instructive to note a relation between the present ap- 
proach and the renormalization group Monte Carlo method 
13(1 Bill used to study phase transitions in lattice models (e.g., 
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polymer models and Ising models for ferromagnets) as well 
as in the quantum field theory. The renormalization procedure 
introduces a scale change in the system. During this procedure 
one consecutively obtains a more and more coarse-grained de- 
scription of the system. Equations and were in fact 
used by Swendsen et al. 13(1 l3lll to describe how the pa- 
rameters of a Hamiltonian change during the renormalization 
procedure. It now appears that the applications of this method 
are more general than the original lattice systems near a phase 
transition, and cover even soft matter problems allowing us to 
"renormalize" Hamiltonians of molecular systems in such a 
way that only the degrees of freedom of primary interest are 
kept in the coarse-grained system. 

C. Dissipative Particle Dynamics (DPD) 

Dissipative particle dynamics was introduced in 1992 by 
Hoogerbrugge and Koelman for simulations of hydrodynamic 
phenomena in complex fluids | lHH. The original formulation 
did not obey detailed balance, though, and in 1995 Espanol 
and Warren |3] formulated a new DPD algorithm which they 
showed to be fully consistent with statistical mechanics. This 
algorithm is nowadays known as DPD. In the following we 
present the basic DPD formalism and discuss some of its fea- 
tures that are relevant to the present work. 

In the basic formulation of DPD, the interactions are pair- 
wise additive and the force exerted by particle j on particle 
i is given as a sum of conservative, dissipative, and random 
forces through py = F^j + F^ + F R , respectively. These 
forces are typically given as 

Ffj = F^ir^Sij, (8) 
Fij = -lu D (rij)(vij ■ ey) ey, (9) 
Pg = <7W fl (r<7)C«3i, (10) 

where r*y = f% — fj is the position vector connecting two par- 
ticles, ry = \fij \ is the interparticle distance, ey = fy/ry 
is the corresponding unit vector, and uy — Vi — Vj is the 
relative velocity of particles i and j. The terms £y are sym- 
metric random variables with zero mean and unit variance, 
and are independent for different pairs of particles and dif- 
ferent times. The conservative force is essentially given by 
Fij {fij)- Constants 7 and a give the amplitudes of the dis- 
sipative and random forces, respectively, and oj d and uj r are 
the corresponding weight functions. 

Let us now return to the different forces. The DPD for- 
mulation does not specify the form of the conservative force. 
The most common choice in standard DPD simulations is 
Fjffa) = ay (1 - ry)/r c for r < r c and Fy c) (ry) = 
otherwise. In other words, the conservative force is derived 
from a soft pair potential U = ay (1 — ry/r c ) 2 /2, where r c 
is the cutoff distance. Interactions between different types of 
coarse-grained particles can then be described by varying the 
amplitude of the conservative interaction ay 1 32] . 

The above is the preferred formulation for the conserva- 
tive force when generic, simple, and soft interactions are ade- 



quate. It was shown by Forrest and Suter in 1995 that coarse- 
graining of molecular representation leads to this type of soft- 
ening 1 33]. However, in many cases this formulation is of too 
generic nature. Although it is possible to use mean field theo- 
ries such as the Flory-Huggins theory for polymers to find the 
amplitudes for the conservative forces between different types 
of particles 1 32], there are cases where this approach does not 
retain enough information about the actual character of the 
different atoms, molecules, or interactions (see discussion be- 
low). 

The weight functions for the dissipative and random forces, 
uj d and uj r , cannot be chosen arbitrarily. This is easy to 
understand by the intuitive argument that the thermal heat 
generated by the random force must be balanced by dissipa- 
tion. Espanol and Warren |3] studied this situation analyti- 
cally and derived a fluctuation-dissipation relation connecting 
the weight functions and amplitudes of the dissipative and ran- 
dom forces via 

u D (r i j) = [u R (r i j)] 2 and a 2 = 2 1 k B T. (11) 

These conditions guarantee convergence to the canonical en- 
semble as required. It is important to notice that oj d and us R 
are completely decoupled from the conservative part. In other 
words, we can consider the DPD algorithm as a momentum 
conserving thermostat for an arbitrary conservative potential. 
This is in fact the way how we apply the DPD algorithm in 
our coarse-grained simulations: The conservative potential is 
defined through the Inverse Monte Carlo procedure and then 
used in the DPD algorithm. 

Finally, to study the dynamics of the system one needs to 
evolve the system in time. In DPD this is simply done by inte- 
grating the Newton's equations of motion. In our simulations 
we used a velocity-Verlet based algorithm adapted for DPD 
(sometimes called DPD-VV) |H0. 

D. Combining MD and DPD 

With the above definitions we now have all the tools for 
coarse-graining. We remove fast internal and orientational de- 
grees of freedom from all water molecules, and represent wa- 
ter molecules as one-site particles interacting with other wa- 
ter particles and ions by spherically symmetric potentials. In 
comparison to the original all-atom model, we have approxi- 
mately one third of the degrees of freedom left. This can be 
interpreted as intermediate level of coarse-graining. 

Without any loss of generality we can choose the weight 
functions to be of the standard form uj R (rij) = 1 — ry/r c 
as discussed above. This choice has been made invariably 
in DPD simulations. Since the Inverse Monte Carlo proce- 
dure provides us with the effective potentials, thus yielding us 
-Fy(r), the only task left is to ensure that the MD and DPD 
systems correspond to each other. In this study, that was done 
by adjusting the amplitude of the dissipative force 7 which de- 
termines the dissipation rate in the DPD system. Here, 7 was 
determined by requiring the short-time decay of the single- 
particle velocity autocorrelation function 

<j>(t) = (*?,(*)• $(0)) (12) 
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0.0 0.1 0.2 0.0 0.1 0.2 0.0 0.1 0.2 
t [ps] t [ps] t [ps] 

FIG. 1: The decay of the single-particle velocity autocorrelation function 4>{t) at early times. Results shown here are for water, Na + and CP 
ions. As the data illustrates, the early-time decay of <f>(t) is essentially identical between MD and DPD simulations for 7 = 0.72 used in the 
DPD simulations. 



to be approximately the same in both MD and DPD systems 
as shown in Fig.^ It shows how the early-time decay of <f>{t) 
can be matched, leading to a value of 7 = 0.72 for the DPD 
model. The long-time behavior of the velocity autocorrelation 
function between MD and DPD is different, however. This is 
an obvious result since some microscopic degrees of freedom 
have been coarse grained and thus effects due to hydrogen 
bonds, for example, are implicitly included in the effective in- 
teractions. This leads to an enhanced diffusion rate at inter- 
mediate times, and is demonstrated in Fig.^as a positive tail 
for the DPD particles. 

This approach is meaningful since (2) it makes sure that 
the early-time dynamics is described properly, while (ii) it 
does not fix the tracer diffusion coefficient. One should no- 
tice that the tracer diffusion coefficient is an integral over the 
velocity correlation function over long times, and in hydro- 
dynamic systems the long-time tail is important. With this 
choice, we can assume the diffusion coefficient to be an inde- 
pendent quantity so that its behavior, found by DPD, can be 
compared with both MD simulations and experiments. Fixing 
7 thus determines the transport properties of the system. 



III. RESULTS 

A. Constructing the potentials 

One of the typical simplifications used in molecular simu- 
lations is the replacement of solvent molecules by continuum 
media. For example, in the primitive electrolyte model, ions 
in water are substituted by charged spheres moving in dielec- 
tric media with the dielectric constant set to about 80. This is 
a serious simplification at small distances (a few angstroms) 
where it is impossible to define a dielectric constant. Besides 
this, in the primitive electrolyte model the ion radius is an ad- 
justable parameter without any clear physical meaning. A bet- 
ter model for effective ion-ion interactions in an aqueous solu- 
tion must take into account the solvation structure around the 
ions. Practically, effective solvent-mediated ion-ion potentials 



may be constructed by the Inverse Monte Carlo method from 
ion-ion radial distribution functions generated in hi gh-q uality 
all-atomic molecular dynamics simulations IU71 12^13411 . This 
is the approach we used in obtaining the effective potentials. 

The MD simulations were performed in the NVT ensem- 
ble using the flexible SPC water model 13511 and Smith-Dang 
parameters for Na + and CI ions 13611 . i.e., the Lennard- 
Jones parameters for the sodium ions were a = 2.35 A and 
e = 0.544 kJ/mol, and for chloride a = 4.4 A and e — 
0.419 kJ/mol. In the results reported here we set the temper- 
ature to 300 K. The electrostatic interactions were taken into 
account by Ewald summation. Other simulation details are 
described in Refs. 1171 13411 . From the MD simulations, the ra- 
dial distribution functions between different pairs of particles 
were calculated and fed as an input in the Inverse Monte Carlo 
procedure. The results are shown in Figs. [2] and [3] 

Figure [2] shows effective potentials between water 
molecules, presented as one-center particles, and between 
other water molecules and ions. For comparison, water-water 
effective potential calculated from MD simulations for pure 
water is also displayed. It is interesting to note that the pres- 
ence of ions has practically no effect on the water-water effec- 
tive potential. 

Figure|3]displays the effective potentials between the ions. 
They are compared to the effective potentials calculated with- 
out water JT^Hl- The effective potentials are very similar 
in the two cases both of them having an oscillating charac- 
ter with one or two oscillations. The potential approaches the 
primitive model potential with dielectric constant close to 80. 
At distances above 10 A the effective potential coincides al- 
most perfectly with the Coulombic potential. 

Within the inverse Monte Carlo procedure, calculation of 
effective potentials without water is much easier than in the 
presence of water. On the other hand, the presence of particles 
representing water is necessary in DPD simulations. Our test 
studies have shown that while the use of effective potentials 
calculated without water may give qualitatively satisfactory 
results, inclusion of water in the Inverse Monte Carlo simula- 
tions essentially leads to an improvement of results, especially 
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FIG. 2: Effective potentials for coarse-grained water in NaCl solu- 
tion. Thin dotted line is water-water effective potential in pure water. 
"O" refers to an oxygen atom in the water molecule. 



B. Coarse-grained simulations 

In our first paper using this approach fiUl . we showed that 
the radial distribution functions computed from the coarse- 
grained DPD simulations for the pairs between the Na + and 
CP ions agreed with those obtained from the all-atom MD 
simulations. In addition to the static properties, the diffusion 
coefficients were shown to be in good agreement between the 
MD and DPD results at various NaCl concentrations between 
0.5 M and 4. 1 M. What makes those results very interesting is 
the fact that, for all the DPD simulations, we used the effec- 
tive potentials obtained at one single NaCl concentration only. 
In that case, the effective potentials were obtained from MD 
simulations at a salt concentration of 0.87 M and then used 
in DPD simulations for various salt concentrations between 
0.5 M and 4.1 M. When compared with MD simulations, the 
agreement was found to be excellent. Here, we present results 
for static correlations, coordination numbers, velocity auto- 
correlation functions, and viscosity. 




FIG. 3: Ion-ion effective potentials in aqueous solution. Solid lines 
describe potentials in presence of one-site water molecules (this 
work), while dashed lines are for effective potential without water 
from Ref. |34|l. 



in the case of dynamical properties. 

It is worth pointing out that the effective potential dif- 
fers from the potential of mean force defined as tppMF — 
—ksT In g(r), which corresponds to the Kirkwood approx- 
imation for the n-particle correlation function. The potential 
of mean force is screened by the other ions in a system and 
decays as {l/r)e^~ r ' rD ' , i.e., not as the 1/r Coulombic po- 
tential. The potential of mean force is therefore not a very 
suitable choice to present ion-ion interactions within a con- 
tinuum solvent model at finite ion concentrations even in a 
coarse-grained solvent system. The effective potentials in- 
clude a contribution from Coulombic forces, however. 





molarity: 


0.55 M 


0.87 M 


2.2 M 


4.1 M 


MD: 


Na+: 


5.8 


5.5 


5.42 


5.3 




CP: 


6.9 


6.8 


6.9 


7.1 


DPD: 


Na+: 


5.67 


5.66 


5.55 


5.40 




CP: 


6.91 


6.96 


6.97 


7.00 



TABLE I: Coordination numbers obtained from MD and DPD sim- 
ulations at various salt concentrations. The error bars for both cases 
are ±0.05. The MD data is from Ref. 1 391. 



It is important to notice that due to the self-consistency 
conditions for the effective potentials, the DPD simulations 
will always produce the correct pair correlation behavior at 
the point in phase space where the effective interactions were 
determined in the first place. As an example, this is demon- 
strated in Fig.|4]for all different kinds of pairs of particles stud- 
ied in this work. It is very relevant to ask, however, how the 
effective interactions change when system variables such as 
the salt concentration or the temperature of the system change. 
We have explored this matter through combined MD and DPD 
simulations and found that for a wide range of salt concentra- 
tions in the present system it is not necessary to recalculate 
the effective potentials. This is demonstrated by the data in 
Table H] which compares the coordination numbers obtained 
from DPD and MD simulations for g(r). In the DPD sim- 
ulations we have again used the effective potentials obtained 
at 0.87 M. We find that both MD and DPD follow the same 
qualitative behavior. For sodium the coordination number de- 
creases with increasing salt concentration whereas for chlo- 
ride we find an opposite trend. Furthermore, the quantitative 
agreement between the MD and DPD results is remarkably 
good. 

Our finding that the effective interactions depend only 
weakly on thermodynamic conditions is also supported by 
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FIG. 4: Radial distribution functions g(r) for different pairs of particles in the aqueous NaCl system studied by both DPD (shown on the left) 
and MD (on the right). The studies were made at a salt concentration of 0.87 M. "O" refers to the ogygen atom in the water molecule. Note 
the agreement between MD and DPD results. 



the present results for shear viscosity (see below) and ear- 
lier studies of the same system for tracer diffusion El. In 
Ref. 0, we found that DPD simulations with effective po- 
tentials yielded tracer diffusion coefficients in good agreement 
with MD simulations. Further support for this result is given 
by previous work |25, 34], where the dependence of effective 
potentials on the salt concentration was studied in detail. It 
was found that effective potentials depend very little on the 
salt concentration: An increase in salt concentration leads to 
a slow decrease in the dielectric permittivity. 

Our results thus suggest that the effective interactions 
change only slightly when the system parameters are varied, 
as is here the case for salt concentration. However, while this 
result holds true in the present system for an aqueous solution 
of monovalent ions, the generality of this finding remains an 
open question. Further studies of divalent systems and more 
complex molecules, among others, are therefore called for. 

Next, we determined the shear viscosity using a Green- 
Kubo relation. As with the other quantities, the shear viscos- 
ity should be compared to results from the MD simulations. 



However, since the shear viscosity coefficient is a collective 
quantity, meaning that all particles in a system give rise to a 
single sample, an accurate determination of the shear viscos- 
ity coefficient through MD simulations is both difficult and 
very time-consuming. Thus it was not done in the present 
case, and as a matter of fact, to the best of our knowledge, 
there are no published reports of MD simulations of shear vis- 
cosity in NaCl solutions. Thus, we compare our DPD data 
to experimental results |37] instead. The results shown in 
Fig.|5]indicate that the qualitative behavior of the shear viscos- 
ity coefficient obtained through DPD simulations is in good 
agreement with the experiments. This result is truly remark- 
able as it shows that both the equilibrium (static) and dynamic 
behavior of the system are properly described by the coarse- 
grained approach. This provides us with strong support that 
the present approach of coarse graining an MD system and ap- 
plying the obtained effective interactions in DPD simulations 
is a promising approach indeed. 

We would like to note that we have very recently performed 
simulations using LiCl and CaCl2 13(811 to study further aspects 
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FIG. 5: The shear viscosity coefficient r\ from coarse-grained 
DPD simulations and experiments. The experimental data is from 
Ref. l37ll . The data by DPD has been scaled to allow comparison. 



of this issue, such as the effects of salt concentration, valence, 
and temperature. Despite the very different natures of these 
solutions (as compared to NaCl), the conclusions with regard 
to the coarse-graining method and its general applicability re- 
main the same. 



IV. CONCLUSIONS 

In this article we have described in detail a method that 
combines the Inverse Monte Carlo method and the momen- 
tum conserving dissipative particle dynamics thermostat for 
coarse-grained simulations. This approach allows systematic 
coarse-graining of molecular systems and conserves hydro- 
dynamics. Here, the level of coarse-graining has been quite 
modest and the potentials used did not have particularly soft 
cores. This is due to the fact that in this system of water and 
NaCl, we coarse-grained the molecular representation of wa- 
ter molecules only. Thus, the corresponding effective interac- 
tions retained features which are typical for systems described 
on an atomic level. However, it is important to notice two 
facts. First of all, further coarse-graining is not possible for 
this kind of a system if one wishes to preserve the identity 
of individual ions and the essential characteristics of water 
molecules. One could, of course, continue the coarse-graining 
procedure to obtain softer potentials by clustering several wa- 
ter molecules together and treating the clusters as new "effec- 
tive" solvent particles. This would lead to further computa- 
tional gains but the ions and molecules would then lose essen- 
tial parts of their identity. As a result of this the interactions 
would become softer and resemble the quadratic potentials 
used for conservative interactions in the original DPD formu- 
lation. Second, despite the modest level of coarse-graining in 



this work, the computational gain is still remarkable, around 
a factor of 20. This speed-up is due to the fact that no explicit 
Ewald summation was required for the coarse-grained system 
(since the effective potentials include information about the 
electrostatics, and the screening length was shorter than the 
applied cut-off) and that SPC water used in molecular dynam- 
ics simulations is replaced by a single coarse-grained particle 
with the electrostatic interactions included in the effective po- 
tentials. 



This approach allows one to coarse-grain complex macro- 
molecular systems to a level where a detailed molecular de- 
scription is replaced with a simplified model of coarse-grained 
particles. For example, a lipid-acyl chain ". . . -CH2 -CH2 - 
CH2-CH2-CH2-CH2 . . ." can be described as ". . .-B- 
B-. . ." where each of the beads "B" represent 2-4 methyl 
groups CH2. Since the beads "B" describe the interaction site 
for a whole molecular group, it is clear that the interactions 
between different beads will be relatively soft in the spirit of 
the original DPD formulation, yet retaining their characteristic 
features well beyond the mean-field description of the original 
DPD formulation. For macromolecular systems, this coarse- 
graining procedure thus has the potential to lead to a major 
reduction in computational effort as well as to softer interac- 
tion potentials. The computational gain depends on the type 
of interactions: While here in the case of NaCl we did not 
need to use Ewald summation, the situation may be different 
in the case of divalent ions such as calcium. In addition, one 
should note that coarse-graining of a macromolecular system 
requires a lot of insight, and there is no unique way of select- 
ing the units to be coarse-grained. Work is in progress to ad- 
dress these issues and the coarse-graining of macromolecular 
systems in general. 



The approach presented in this article has close connections 
to other methods as well. The Inverse Monte Carlo procedure 
has the same foundations as the integral equation methods of 
Bolhuis et al. I KIM ill . Second, the coarse-graining procedure 
used by Reith, Mayer and Miiller-Plathe 11511 has a close re- 
semblance to the Inverse Monte Carlo procedure. The latter 
group has also applied their method to coarse-grain polymeric 
systems in the spirit described in the previous paragraph. The 
novelty here is the application of the effective potentials in 
coarse-grained simulations by using the dissipative particle 
dynamics thermostat. The results for the present system in- 
dicate that dynamic properties could be studied using this ap- 
proach. Whether or not this holds true in more complex sys- 
tems, and what are the limits of this approach are topics of 
further studies. 



In conclusion, the agreement between the DPD simula- 
tions, the MD results and experiments strongly indicates that 
the method presented in this work is both physically sound 
and provides a controllable method for performing coarse- 
graining and mesoscale simulations at least to a moderate de- 
gree. 
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